perm filename COMFM.F4[JC,MUS]1 blob sn#074713 filedate 1973-11-26 generic text, type T, neo UTF8
00100	C This is taken from the "Handbook of Mathematical Functions"
00200	C by Abramowitz and Stegun, Dover press S1272, 1965.
00300	C Chapter 9 is on Bessel Functions. The polynomial
00400	C approximations are found in section 9.4, pages 369
00500	C and 370. The recursion formula for generating higher
00600	C orders from J0 and J1 is found at the bottom of
00700	C table 9.1, page 390
00800	
00900		SUBROUTINE BESCOEF(MI,X,I)
01000		REAL MI,J0,J1
01100		COMMON FREQ(3,0/200)
01200		IF(I.GT.1)GO TO 30
01300		IF(MI.GE.3.0)GO TO 10
01400		xs=(MI/3)**2
01500		j0=1.0+xs*(-2.2499997+xs*(1.2656208+xs*(-0.3163866+
01600		1xs*(0.0444479+xs*(-0.0039444+xs*0.00021)))))
01700		j1=MI*(.5+xs*(-0.56249985+xs*(0.21093573+
01800		1xs*(-0.03954289+xs*(0.00443319+xs*(-0.00031761+
01900		2xs*0.00001109))))))
02000		GO TO 20
02100	10	xs=3.0/MI
02200		j0=(1.0/sqrt(MI))*(0.79788456+xs*(-0.00000077+
02300		1xs*(-0.0055274+xs*(-0.00009512+xs*(0.00137237+
02400		2xs*(-0.00072805+xs*0.00014476))))))*
02500		3cos(MI-0.78539816+xs*(-0.04166397+
02600		4xs*(-0.00003954+xs*(0.00262573+
02700		5xs*(-0.00054125+xs*(-0.00029333+
02800		6xs*0.00013558))))))
02900		j1=(1.0/sqrt(MI))*(0.79788456+xs*(0.00000156+
03000		1xs*(0.01659667+
03100		1xs*(0.00017105+xs*(-0.00249511+xs*(0.00113653-
03200		2xs*0.00020033))))))*
03300	 	3cos(MI-2.35619449+xs*(0.12499612+xs*(0.0000565+
03400		4xs*(-0.00637879+xs*(0.00074348+xs*(0.00079824-
03500		5xs*0.00029166))))))
03600	20	FREQ(2,0)=J0
03700		FREQ(2,1)=J1
03800		FREQ(2,2)=J1
03900		RETURN
04000	30	Y=I
04100		IF(MI.LE.0.0000001)GO TO 50
04200		IJ=I*2		
04300		X=((2.0*(Y-1.0))/MI)*FREQ(2,IJ-2)-FREQ(2,IJ-4)
04400		RETURN
04500	50	X=0
04600		RETURN
04700		END
     

00100		COMMON FREQ(3,0/200)
00150		IT=-1.0
00200		ACCEPT 102,CAR,XMOD,CINDEX
00300		MAX=CINDEX+5.0
00400		DO 100 J=0,MAX*2,2
00500		I=J/2
00600		NORDER=I
00700		IF(I.NE.1)CALL BESCOEF(CINDEX,COEF,NORDER)
00800		IF(I.GT.0)GO TO 20
00900		FREQ(1,I)=CAR
01000		FREQ(3,I)=0.0
01100		GO TO 30
01200	20	XI=I
01300		YMOD=XI*XMOD
01400		FREQ(1,J-1)=CAR+YMOD
01500		FREQ(1,J)=CAR-YMOD
01600		IF(I.EQ.1)GO TO 10
01700		FREQ(2,J-1)=COEF
01800		FREQ(2,J)=COEF
01900	10	FREQ(3,J-1)=I
01950		IF(IT)I=-I
02000		FREQ(3,J)=I
02050		IT=-IT
02100	30	CONTINUE
02200	100	CONTINUE
02300		TYPE 101,((FREQ(K,L),K=1,3),L=0,MAX*2)
02400	101	FORMAT(1H 20(3F/))
02500	102	FORMAT(3F)
02600		END